INCLUDE 'mkl_dss.f90'
module MKL_DSS_Solver_m
use SMatrix_m
use mkl_dss
implicit none
INTEGER, PARAMETER :: dp = KIND(1.0D0)
type MKL_DSS_Solver
    logical :: isSym
    integer :: error
    integer :: nRhs
    integer :: nRows
    integer :: nCols
    integer :: nNonZeros
    REAL(KIND=DP), ALLOCATABLE :: values( : )
    integer, allocatable :: rowIndex(:), columns(:)
    REAL(KIND=DP), ALLOCATABLE :: statOUt( : )
    CHARACTER*15 :: statIn
    integer :: perm(1)
    integer :: buff(20)
    type(MKL_DSS_HANDLE) :: handle
contains
    procedure :: MKLDSS_initSolver
    procedure :: MKLDSS_initSolver2
    generic :: initSolver => MKLDSS_initSolver, MKLDSS_initSolver2
    procedure :: MKLDSS_Solve
    procedure :: MKLDSS_SolveOne
    generic :: solve => MKLDSS_Solve, MKLDSS_SolveOne
    procedure :: printStatistics => MKLDSS_printStatistics
    procedure :: printErrorInfo => MKLDSS_printErrorInfo
    procedure :: disposeSolver => MKLDSS_disposeSolver
end type
private :: MKLDSS_initSolver, MKLDSS_Solve, MKLDSS_printStatistics, MKLDSS_printErrorInfo, MKLDSS_disposeSolver
contains
    
subroutine MKLDSS_initSolver(this, A)
    implicit none
    class(MKL_DSS_Solver), intent(out) :: this
    type(SMatrix), intent(inout) :: A
    this.nRows = A%M
    this.nCols = A%N
    this.perm(1) = 0
    if ( .NOT. A.isSorted) call A.sort()
    this.isSym = A.isSymmetric()
    call A.exportCSRFormat(this.isSym, this.rowIndex, this.columns, this.values)
    this.nNonZeros = size(this.columns, 1)
    this.error = DSS_CREATE(this.handle, MKL_DSS_DEFAULTS)
    IF (this.error /= MKL_DSS_SUCCESS)  call this.printErrorInfo()
    if (this.isSym) then
        this.error = DSS_DEFINE_STRUCTURE(this.handle, MKL_DSS_SYMMETRIC, this.rowIndex, A%M, &
                                 A%N, this.columns, this.nNonZeros)
    else
        this.error = DSS_DEFINE_STRUCTURE(this.handle, MKL_DSS_NON_SYMMETRIC, this.rowIndex, A%M, &
                                 A%N, this.columns, this.nNonZeros)
    end if
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()

    write(*, *) "Reorder the matrix...."
    this.error = DSS_REORDER(this.handle, MKL_DSS_DEFAULTS, this.perm)
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
    write(*, *) "Factor the matrix...."
    this.error = DSS_FACTOR_REAL(this.handle, MKL_DSS_DEFAULTS, this.values)
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
end subroutine

subroutine MKLDSS_initSolver2(this, A, isSym, isDefinite)
    implicit none
    class(MKL_DSS_Solver), intent(out) :: this
    type(SMatrix), intent(inout) :: A
    logical, intent(in) :: isSym, isDefinite
    
    this.nRows = A%M
    this.nCols = A%N
    this.perm(1) = 0
    if ( .NOT. A.isSorted) call A.sort()
    call A.exportCSRFormat(isSym, this.rowIndex, this.columns, this.values)
    this.nNonZeros = size(this.columns, 1)
    this.error = DSS_CREATE(this.handle, MKL_DSS_DEFAULTS)
    IF (this.error /= MKL_DSS_SUCCESS)  call this.printErrorInfo()
    if (isSym) then
        this.error = DSS_DEFINE_STRUCTURE(this.handle, MKL_DSS_SYMMETRIC, this.rowIndex, A%M, &
                                 A%N, this.columns, this.nNonZeros)
    else
        this.error = DSS_DEFINE_STRUCTURE(this.handle, MKL_DSS_NON_SYMMETRIC, this.rowIndex, A%M, &
                                 A%N, this.columns, this.nNonZeros)
    end if
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()

    write(*, *) "Reorder the matrix...."
    this.error = DSS_REORDER(this.handle, MKL_DSS_DEFAULTS, this.perm)
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
    write(*, *) "Factor the matrix...."
    if (isDefinite) then
        this.error = DSS_FACTOR_REAL(this.handle, MKL_DSS_DEFAULTS, this.values)
    else
        this.error = DSS_FACTOR_REAL(this.handle, MKL_DSS_INDEFINITE, this.values)
    end if
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
end subroutine

subroutine MKLDSS_Solve(this, b, n, x, verbosity)
    implicit none
    class(MKL_DSS_Solver) :: this
    real*8, intent(inout) :: b(:), x(:)
    integer, intent(in) :: n, verbosity
    
    integer :: i
    this.error = DSS_SOLVE_REAL(this.handle, MKL_DSS_DEFAULTS, b, n, x)
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
    if (verbosity >= 1) then
        do i = 1, n
            write(*, *) i, "min", minval(x((i-1)*this.nRows+1:i*this.nRows)), &
                           "max", maxval(x((i-1)*this.nRows+1:i*this.nRows))   
        end do
    end if
    if (verbosity >= 2) call this.printStatistics()
end subroutine

subroutine MKLDSS_SolveOne(this, b, x, verbosity)
    implicit none
    class(MKL_DSS_Solver) :: this
    real*8, intent(inout) :: b(:), x(:)
    integer, intent(in) :: verbosity
    
    this.error = DSS_SOLVE_REAL(this.handle, MKL_DSS_DEFAULTS, b, 1, x)
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
    if (verbosity >= 1) write(*, *) "min", minval(x), "max", maxval(x)
    if (verbosity >= 2) call this.printStatistics()
end subroutine

subroutine MKLDSS_disposeSolver(this)
    implicit none
    class(MKL_DSS_Solver) :: this
    this.error = DSS_DELETE(this.handle, MKL_DSS_DEFAULTS)
    IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
    IF (ALLOCATED(this.rowIndex)) DEALLOCATE(this.rowIndex)
    IF (ALLOCATED(this.columns )) DEALLOCATE(this.columns)
    IF (ALLOCATED(this.values )) DEALLOCATE(this.values)
    IF (allocated(this.statout)) deallocate(this.statout)
end subroutine

subroutine MKLDSS_printStatistics(this)
    implicit none
    class(MKL_DSS_Solver), intent(inout) :: this
    IF(this.nRows < this.nNonZeros) THEN
       if (.NOT.allocated(this.statOut)) ALLOCATE(this.statOut(5))
       this.statIn = 'determinant'
       call mkl_cvt_to_null_terminated_str(this.buff,20,this.statIn)
       this.error = DSS_STATISTICS(this.handle, MKL_DSS_DEFAULTS, this.buff, this.statOut)
       IF (this.error /= MKL_DSS_SUCCESS) call this.printErrorInfo()
       WRITE(*,"('pow of determinant is '(5F10.3))") (this.statOut(1))
       WRITE(*,"('base of determinant is '(5F10.3))") (this.statOut(2))
       WRITE(*,"('Determinant is '(5F10.3))") ((10**this.statOut(1))*this.statOut(2))
    END IF
end subroutine

subroutine MKLDSS_printErrorInfo(this)
    implicit none
    class(MKL_DSS_Solver), intent(in) :: this
    WRITE(*,*) "Solver returned error code ", this.error
    STOP 1    
end subroutine

end module